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ABSTRACT 

We introduce VErsatile SPectral Analysis (VESPA): a new method which aims to 
recover robust star formation and metallicity histories from galactic spectra. VESPA 
uses the full spectral range to construct a galaxy history from synthetic models. We 
investigate the use of an adaptative parametrization grid to recover reliable star for- 
mation histories on a galaxy- by-galaxy basis. Our goal is robustness as opposed to high 
resolution histories, and the method is designed to return high time resolution only 
where the data demand it. In this paper we detail the method and we present our find- 
ings when we apply VESPA to synthetic and real Sloan Digital Sky Survey (SDSS) 
spectroscopic data. We show that the number of parameters that can be recovered 
from a spectrum depends strongly on the signal-to-noise, wavelength coverage and 
presence or absence of a young population. For a typical SDSS sample of galaxies, we 
can normally recover between 2 to 5 stellar populations. We find very good agreement 
between VESPA and our previous analysis of the SDSS sample with MOPED. 

Key words: methods: data analysis - methods: statistical - galaxies: stellar content 
- galaxies: evolution - galaxies: formation 



1 INTRODUCTION 

The spectrum of a galaxy holds vasts amounts of informa- 
tion about that galaxy's history and evolution. Finding a 
way to tap directly into this source of knowledge would 
not only provide us with crucial information about that 
galaxy's evolutionary path, but would also allow us to 
integrate this knowledge over a large number of galaxies 
and therefore derive cosmological information. 

Galaxy formation and evolution are still far from being 
well understood. Galaxies are extremely complex objects, 
formed via complicated non-linear processes, and any 
approach (be it observational, semi-analytical or compu- 
tational) inevitably relies on simplifications. If we try to 
analyse a galaxy's luminous output in terms of a history 
parametrized by some chosen physical quantities, such a 
simplification is also in order. The reason is two-fold: firstly 
we are limited by our knowledge and ability to model 
all the physical processes which happen in a galaxy and 
produce the observed spectrum we are analysing; secondly, 
the observed spectrum is inevitably perturbed by noise, 
which intrinsically limits the amount of information we can 
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Measuring and understanding the star formation history 
of the Universe is therefore essential to our understanding 
of galaxy evolution - when, where and in what conditions 
did stars form throughout cosmic history? The traditional 
and simplest way to probe this is to measure the observed 
instantaneous star formation rate in galaxies at differ- 
ent redshifts. This can be achieved by looking at light 
emitted by young stars in the ultra-violet (UV) band or 
its s e condary effects , (e.g . iMadau et al.l [l996; Kcn nicutt 
1 19981 : iHopkins et all l200fj: iBundv et all 120061; lErb et at 
huod: I Abraham et alj|2007t iNoeske et al.ll2007l; IVerma et al 



l2007f ). A complementary method is to look at present day 
galaxies and extract their star formation history, which 
spans the lifetime of the galaxy. Different teams have 
analysed a large number of galaxies in this way, whether 



20031: Panter et all 20031; Cid Fernandes et all 


2004; 


Heavens et alj|2004 iMathis et al.ll2006l; 


Ocvirk et al.l 


2006; 


Panter et al.l 120061: ICid Fernandes et al. 


l2007t). or by con- 



centrating on particular spectral features or indices (e.g. 



Kauff mann et al. 120031 : Tremonti et al. 2004 iGallazzi et al.l 
20051 : iBarber et al. 2006), which are known to be correlated 



with age or metallicities (e.g. IWorthevlll994l : iThomas et al.l 
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l2003h . 

To do this, we rely on synthetic stellar population models to 
describe a galaxy in terms of its stellar components, but by 
modelling a galaxy in this way we are intrinsically limited 
by the quality of the models. There are also potential 
concerns with flux calibration errors. However, using the 
full spectrum to recover the fossil record of a galaxy - or of 
an ensemble of galaxies - is an extremely powerful method, 
as the quality and amount of data relating to local galaxies 
vastly outshines that which concerns high-redshift galaxies. 
Splitting a galaxy into simple stellar populations of different 
ages and metallicities is a natural way of parameterising 
a galaxy, and it allows realistic fits to real galaxies (e.g. 
iBruzual fc Charlotll2003l ). Galactic archeology has become 
increasingly popular in the literature recently, largely due to 
the increase in sophistication of stellar population synthesis 
codes and the improvement of the stellar spectra libraries 
upon which they are based, and also due to the availability 
of large well-calibrated spectroscopic d atabases, such as 
the Sloan Digital Sky Survey (SDSS) l|York et all |2000| ; 
IStrauss et afl|2002h . 

In any case, without imposing any constraints on the 
allowed form of the star formation history, or perhaps an 
age-metallicity relation, the parameter space can become 
unsustainably large for a traditional approach. Ideally, one 
would like to do without such pre-constraints. Recently, 
different research teams have come up with widely different 
solutions for this pr oblem. MOPED Jlleavens et aJj|2000h 
and STARLIGHT (|Cid Fernandes et all 120041 ) explore a 
well-chosen parameter space in order to find the best 
possible fit to the data. In the case of MOPED, this relies 
on compression of the full spectrum to a much smaller set 
of numbers which retains all the information about the 
parameters it tries to recover; STARLIGHT on the other 
hand, searches for its best fit using th e full spectrum with 
a Metropolis algorithm. STECMAP (jOcvirk et al.l 120061 ) 
solves the problem using an algebraic least-squares solution 
and a well-chosen regularization to keep the inversions 
stable. All of these and other methods acknowledge the 
same limitation - noise in the data and in the models 
introduces degeneracies into the problem which can lead 
to unphysical results. MOPED, for example, has produced 
some remarkable results concerning the average star forma- 
tion history of the Universe by analysing a large sample 
of galaxies. However, MOPED's authors have cautioned 
against over-interpreting the results on a galaxy-by-galaxy 
basis, due to the problem mentioned above. This is directly 
related to the question of how finely one should param- 
eterise a galaxy, and what the consequences of this might be. 

Much of the motivation for VESPA came from the reali- 
sation that this problem will vary from galaxy to galaxy, 
and that the method of choosing a single parametrization 
to analyse a large number of galaxies can be improved on. 

VESPA is based on three main ideas, which we present here 
and develop further in the main text: 



from any given set of data, and the amount of information 
which can be recovered from an individual galaxy varies. 

• The recovered star formation fractions should be posi- 
tive. 

• Even though the full unconstrained problem is non- 
linear, it is piecewise linear in well-chosen regions of pa- 
rameter space. 

VESPA's ultimate goal is to derive robust information 
for each galaxy individually, by adapting the number of 
parameters it recovers on a galaxy- by- galaxy basis and 
increasing the resolution in parameter space only where the 
data warrant it. In a nutshell, this is how VESPA works: we 
estimate how many parameters we can recover from a given 
spectrum, given its noise, shape, spectral r esolution and 
wavele ngth range using an analysis given by lOcvirk et al.l 
(2006). In that paper, Singular Value Decomposition (SVD) 
is used to find a least squares solution, and this solution 
is analysed in terms of its singular vectors. VESPA uses 
this method only as an analysis of the sol ution, and uses 
Boun ded- Variable Least-Squares (BVLS) (IStark fc Parker! 
Il995l) to reach a non-negative solution in several regimes 
where linearity applies. 

This paper is organised as follows: in Section [2] we present 
the method, in Section [3] we apply VESPA to a variety of 
synthetic spectra, in SectionUwe apply VESPA to a sample 
of galaxies from the Sloan Digital Sky Survey spectroscopic 
database and we compare our results to those obtained with 
MOPED, and finally in Section[S]we present our conclusions. 



2 METHOD 

In this section we lay down the problem to solve in detail, 
and explain the different steps VESPA uses to find a solution 
for each galaxy. 



2.1 The problem 

We assume a galaxy is composed of a series of simple stel- 
lar populations (SSP) of varying ages and metallicities. The 
unobscured rest frame luminosity per unit wavelength of a 
galaxy can then be written as 



Fx = 



1>(t)S*(t,Z)dt 



(1) 



where ip(t) is the star formation rate (solar masses formed 
per unit of time) and S\(t,Z) is the luminosity per unit 
wavelength of a single stellar population of age t and metal- 
licity Z, per unit mass. The dependency of the metallicity on 
age is unconstrained, turning this into a non-linear problem. 

In order to solve this problem, we start by discretizing in 
wavelength and time, by averaging these two quantities 
into well chosen bins. For now we present the problem with 
a generalised parametrization, and discuss our choice in 
Section \2. 31 We will use greek indices to indicate time bins, 
and roman indices to indicate wavelength bins. 



• There is only so much information one can safely recover 



The problem becomes 



Recovering galaxy histories using VESPA 3 



Fj = Vi G(Z ) 



(2) 



where Fj — (Fi, ...,Fd) is the luminosity of the jth wave- 
length bin of width AA, G(Z a ) a j is the jth luminosity point 
of a stellar population of age t a = (ii, ...,ts) (spanning an 
age range of At a ) and metallicity Z a , and x a — (xi, xs) 
is the total mass of population G(Z) a j in the time bin At a . 

Although the full metallicity problem is non-linear, interpo- 
lating between tabulated values of Z gives a piecewise linear 
behaviour: 

(3) 

and the problem then becomes 

Fj = Va; Q [g a G(Z a ] (4) 



where G(Z a ^ a ) a j and G(Zb, a )aj are equivalent to G(Z a ) a j 
as above, but at fixed metallicities Z a and Zb, which bound 
the true Z. If this interpolation matches the models' reso- 
lution in Z, then we are not degrading the models in any way. 

Solving the problem then requires finding the correct metal- 
licity range. One should not underestimate the complexity 
this implies - trying all possible combination of consecutive 
values of Z a and Zb in a grid of 16 age bins would lead to a 
total number of calculations of the order of 10 9 , which is un- 
feasible even in today's fast personal workstations. We work 
around this problem using an iterative approach, which we 
describe in Section \2. 3. 21 



2.1.1 Dust extinction 

An important component when describing the luminous out- 
put of a galaxy is dust, as different wavelengths are affected 
in different ways. The simplest possible approach is to use 
one-parameter dust model, according to which we apply 
a single dust screen to the combined luminosity of all the 
galactic components. Equation (fJJ becomes 

Fx = f dust (r x ) [ ip(t)S x (t,Z)dt (5) 
Jo 

where we are assuming the dust extinction is the same for 
all stars, and characterised by the optical depth, t\. 

However, it is also well known that very young stars are 
likely to be more affected by dust. In an attempt to in- 
clude this in our modelling, we foll ow the two-parameter 
dust model of ICharlot fc Falll (|2000h in which young stars 
are embebbed in their birth cloud up to a time tsci when 
they break free into the inter-stellar medium (ISM): 



Fx = I f duat {Tx,t)i>{t)S x {t,Z)dt 
Jo 

and 

fdust(T\ SM )fdust(T\ C ),t < tBC 
fdust{r{ SM ),t > tBC 



fdust{T\,t) — <. /ISM\ 



(6) 



(7) 



where t{ sm is the optical depth of the ISM and rf c is the 
optical depth of the birth cloud. Following ICharlot fc Falll 



l|2000l ). we take t B c = 0.03 Gyrs. 

There is a variety of choices for the form of fdustijx). To 
model the dust in the ISM, we use the mixed slab model of 
ICharlot fc Fan] l|2000l ) for low optical depths (tv < 1), for 
which 



fdv.st(rx) = - — [1 + (ta 
2ta 



l)eXp(-TA)-Tf£i(TA)] 



(8) 



where E\ is the exponential integral and ta is the optical 
depth of the slab. This model is known to be less accurate 
for high dust values, and for optical depths greater than one 
we take a uniform screening model with 



fdust(rx) = exp(-T A ). 



(9) 



We only use the uniform screening model to model the dust 
in the birth cloud and we use ta = tv(A/5500A) -0 ' 7 as our 
extinction curve for both environments. 

As described, dust is a non-linear problem. In practice, we 
solve the linear problem described by equation Q with a 
number of dust extinctions applied to the matrices G(Z)ij 
and choose the values of t v sm and Ty ° which result in the 
best fit to the data. 



We initially use a binary chop search for r v SM £ [0, 4] and 
keep Ty° fixed and equal to zero, which results in trying out 
typically around nine values of t v sm . If this initial solution 
reveals star formation at a time less than tsc we repeat 
our search on a two-dimensional grid, and fit for Ty SM and 
Ty C simultaneously. There is no penalty except in CPU time 
to apply the two-parameter search, but we find that this 
procedure is robust (see section 3.4). 



2.2 The solution 

In this section we describe the method used to reach a 
solution for a galaxy, given a set of models and a gener- 
alised parametrization. The construction of these models 
and choice of parameters is described in Sections l2.3l and l2.4l 

We re-write the problem described by equation in a sim- 
pler way 



Fj=^c K A Kj {Z K 



(10) 



where Z K — Z a for k < S and Z K = Zt for k > S. A 
is a D x 2S matrix composed of synthetic models at the 
corresponding metallicities, and c = (ci, C2s) is the solu- 
tion vector, from which the x a and g a in equation @ can 
be calculated. We can then calculate a linearly interpolated 
metallicity at age t a 



Z a = g a Z a + (1 - g a )Z b . 



(11) 



For every age t a we aim to recover two parameters: x a - 
the total mass formed at that age (within a period At a ) - 
and Z a - a mass-weighted metallicity. 

At this stage we are not concerned with our choice of t a and 
At a - although these are crucial and will be discussed later. 
For a given set of chosen parameters, we find c, such that 
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x 2 = (£>-£„ c «Ay 



(12) 



is minimised (where <7j is the error in the measured flux bin 
Fj). 

A linear problem with a least squares constraint has a simple 
analytic solution which, for constant oj (white-noise) is 



CLS 



(A T - A)- 1 -A T -F 



(13) 



In principle, any matrix inversion method, e.g. Singular 
Value Decomposition (S VD) , can be used to solve JT3} ■ How- 
ever, we would like to impose positivity constraints on the 
recovered solutions. Negative solutions are unphysical, but 
unfortunately common in a problem perturbed by noise. 



2.2.2 Noise 

The inversion in equation (|13|l is often highly sensitive to 
noise, and care is needed when recovering solutions with 
matrix inversion methods. The fit in data-space will always 
improve as we increase the number of parameters, but these 
might not all pr ovide meaningful in formation. We follow an 
analysis given in lOcvirk et al.l ( |2006h in order to understand 
how much this affects our results, and to choose a suitable 
age parametrization for each galaxy. This is not an exact 
method, and it does not guarantee that the solutions we 
recover have no contribution from noise. However, we found 
that in most cases it provides a very useful guideline (see 
section [3731 in particular Figure [TT}. 

We refer the reader to the above paper for a full discussion, 
and we reproduce here the steps used in our analysis. 



2.2.1 BVLS and positivity 

We use Bounded- Variable Least-Squares (BVLS) 
l|Stark fe Parker! Il995l ) in order to solve CCH- BVLS is 
an algorithm which solves linear problems whose variables 
are subject to upper and lower constraints. It uses an 
active set strategy to choose sets of free variables, and 
uses QR matrix decomposition to solve the unconstrained 
least-squares problem of each candidate set of free variables 
using (|13p: 



cls = (E T •E)- 1 -E T -F 



(14) 



where E is effectively composed of those columns of A 
for which ct is unconstrained, and of zero vectors for 
those columns for which is set to zero. BVLS is an 
extension of the Non-N egative Least Squares algorithm 
|Lawson fc Hanson] 1 19741 ). and they are both proven to 
converge in a finite number of iterations. Positivity is the 
only constraint in VESPA's solution. 

BVLS and positivity have various advantages. Most obvious 
is the fact that we do away with negative solutions. In a 
non-constrained method (such as SVD) negative values are 
a response to the fact that the data is noisy. Similarly, we 
find that zero values returned by BVLS (in, for example, 
a synthetic galaxy with continuous star formation across 
all time) are also an artifact from noisy data. It should be 
kept in mind that, if the method is unbiased, this problem 
is solved by analysing a number of noisy realisations of the 
original problem - what we find is that the true values of 
the parameters we try to recover are consistent with the 
distributions yielded by this process. In this sense, not even 
a negative value presents a problem necessarily, as long as 
it is consistent with zero (or the correct solution). Given 
that we have found no bias when using BVLS, we feel it is 
an advantage to discard a priori solutions we know to be 
unphysical. 

Another advantage to using BVLS is the fact that, by fixing 
some parameters to the lower boundary (zero, in this case), 
it effectively reduces the number of fitting parameters to 
the number of those which keeps unconstrained. Given the 
overall aims of VESPA, this has proven to be advantageous. 



We use SVD to decompose the model matrix E as 
E = U W V T 



(15) 



where U is a D x 25 orthonormal matrix with singular data- 
vectors u K as columns, V is a 2S x 2S orthonormal matrix 
with the singular solution-vectors v K as columns, and W is 
a 2S x 2S diagonal matrix W = diag(w\, W2s) where w K 
are the matrix singular values in decreasing order. Replacing 
E by this decomposition in equation (|13|l gives 



c LS = V- W _1 -U T -F = 



^ W K 



(16) 



The solution vector is a linear combination of the solution 
singular values, parametrized by the dot product between 
the data and the corresponding data singular vector, and 
divided by the k th singular value. The data vector itself is a 
combination of the true underlying emitted flux and noise: 
F = Ftme + e - Equation (|16[) becomes 



2S 



25 



Eu K -F true \ u K -e 
v K + > v« 
Wk, z — ' W K 

K=l K=l 



(17) 



where ctrue is the solution vector to the noiseless problem 
and c e is an unavoidable added term due to the presence of 
noise. 

It is extremely informative to compare the amplitudes 
of the two terms in the sum 1)170 . and to monitor their 
contributions to the solution vector with varying rank. In 
Figure [T] we plot |uj ■ F| and ujfe as a function of rank 
k, for a synthetic spectrum with a SNR per pixel of 50 
(at a resolution of 3A) and an exponentially-decaying star 
formation history . We observe the behaviour described 
and discussed in lOcvirk et alj ([2006). The combinations 
associated with the noise terms maintain a roughly constant 
power across all ranks, with a an average value of (F) /SNR. 
The data terms, however, drop significantly with rank, and 
we can therefore identify two ranges: a noise-dominated 
K-range, in which the noise contributions match or dom- 
inate the true data contributions, and a data-dominated 
range, where the contributions to the solution are largely 
data motivated. We call the transition rank n cr it. Overall, 
high-K ranks tend to dominate the solution, since the 
singular values w K decrease with k. This only amplifies 
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Figure 1. The behaviour of the singular values with matrix rank 
k. The stars are |u^-F| and the squares are u^-e. The line is 
(F) /SNR, which in this case has a value of approximately 0.06. 



the problem by giving greater weight to noise-dominated 
terms in the sum (|16[) . Figure [2] shows the contribution 
coming from each rank k to the final solution - the co- 
efficient (ul-F)/w K . We see this weight increases with rank. 

Whereas this analysis gives us great insight into the 
problem, we do not in fact use the sum (|16[1 to obtain c^s, 
for the reasons given in section [2. 2. II 

For real data we are only able to calculate u^-F and estimate 
the noise level at (F) /SNR and we use this information 
to estimate the number of non-zero parameters to recover 
from the data. Our aim is to a have a solution which is 
dominated by the signal, and not by the noise. We therefore 
want our number of non-zero recovered parameters to be 
less than or equal to K cr it- Estimating where this transition 
happens is always a noisy process. In this paper we take 
the conservative approach of setting n cr it to be the rank at 
which the perturbed singular values first cross the (F) /SNR 
barrier. In the case of Figure [T] this happens at K cr u = 7. 



100 








^ -100- 



-200 - 



-300 r 

5 10 15 20 25 

Rank k 

Figure 2. The coefficients in sum II16I I as a function of rank k. 
We see that the highest rank modes (corresponding to the smaller 
singular values) tend to contribute the most to the solution. 



up to 14 Gyr. The grid has three further resolution levels, 
where we split the age of the Universe in eight, four and 
finally two age bins, also logarithmically spaced in the same 
range. 

The idea behind the multi-resolution grid is to start our 
search with a low number of parameters (in coarser resolu- 
tion, so that the entirety of the age of Universe is covered), 
and then increase the resolution only where the data war- 
rant it by splitting the bin with the highest flux contribution 
in two, and so on. In effect, we construct one such grid for 
each of the tabulated metallicities, Z a and Zi,. We work with 
five metallicity values, Z = [0.0004, 0.004, 0.008, 0.02, 0.05] 
which correspond to the metallicity resolution of the models 
used, where Z is the fraction of the mass of the star 
composed of metals (Zq = 0.02). The construction of the 
models for each of the time bins is discussed in Section 12.41 

To each of the grids we can apply a dust extinction as ex- 
plained in Section \2. 1.1 1 



2.3 Choosing a galaxy parametrization 

One of the advantages of VESPA is that it has the ability 
to choose the number of parameters to recover in any 
given galaxy. This is possible due to a time grid of varying 
resolutions, which VESPA can explore to find a solution. 
This section describes this grid and the criteria used to 
reach a final parametrization. 

2.3.1 The grid 

We work on a grid with a maximum resolution of 16 age 
bins, logarithmically spaced in lookback time from 0.02 



2.3.2 The search 

We go through the following steps in order to reach a solu- 
tion: 

(i) We begin our search with three bins: two bins of width 
4 and one bin of width 8 (oldest), where here we are mea- 
suring widths in units of high-resolution bins. 

(ii) We calculate a solution using equation (|10p for every 
possible combination of consecutive boundaries Z a and Zb, 
and we choose the one which gives the best value of reduced 

x 2 . 

(iii) We calculate the number of perturbed singular values 
above the noise level, as described at the end of Section ^. 2, 21 
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(iv) We find the bin which contributes the most to the 
total flux and we split it into two. 

(v) We find a solution in the new parametrization, this 
time by trying out all possible combinations of Z a and Z\, for 
the newly split bins only, and fixing the metallicity bound- 
aries of the remainder bins to the boundaries obtained in 
the previous solution. If a bin had no stars in the previous 
iteration, we set Z a = 0.0004 and Zb = 0.05. 

(vi) We return to (iii) and we proceed until we have 
reached the maximum resolution in all populated bins. 

(vii) We look backwards in our sequence of solutions for 
the last instance with a number of non-zero recovered pa- 
rameters equal to or less than K C rit as calculated in (iii) and 
take this as our best solution. 

We illustrate this sequence in Figure [3] where we show the 
evolution of the search in a synthetic galaxy composed of 
two stellar bursts of equal star formation rates - one young 
and one old. VESPA first splits the components which con- 
tribute the most to the total flux. In this case this is the 
young burst which can be seen in the first bin. Even though 
VESPA always resolves bins with any mass to the possible 
highest resolution, it then searches for the latest solution 
which has passed the SVD criterion explained in Section 
12.2.21 In this case, this corresponds to the fifth from the top 
solution. VESPA chooses this solution in favour of the fol- 
lowing ones due to the number of perturbed singular values 
above the solid line (right panel). In this case, the solution 
chosen by VESPA is a better fit in parameters space (note 
the logarithmic scale in the y-axis - the following solution 
put the vast majority of the mass in the wrong bin). We ob- 
served this type of improvement in the majority of all cases 
studied (see Figure [TTJ. 



2.3.3 The final solution 

Our final solution comes in a parametrization such that the 
total number of non-zero recovered parameters is less than 
or equal to the number of perturbed singular values above 
the estimated noise level. 

The above sequence is performed for each of several combi- 
nations of Ty C ,Ty SM , and we choose the attenuation which 
provides the best fit. 



2-4-1 High-resolution age bins 

At our highest resolution we work with 16 age bins, equally 
spaced in a logarithmic time scale between now and the 
age of the Universe. In each bin, we assume a constant star 
formation rate 



f™(\,Z)=Tp S(\,t,Z)dt 

JAt a 

with ip = 1/At a . 



2.4-2 Low-resolution age bins 



(18) 



As described in Section \2. 3. II we work on a grid of different 
resolution time bins and we construct the low resolution bins 
using the high resolution bins described in Section ^. 4. II We 
do not assume a constant star formation rate in this case, as 
in wider bins the light from the younger components would 
largely dominate over the contribution from the older ones. 
Instead, we use a decaying star formation history, such that 
the light contributions from all the components are compa- 
rable. Recall equation Q 



fZ R (\,z) = 



ip(t)S(\,t,Z)dt, 



which we approximate to 



Y, a ^fS R {\z)^ a At a 



(19) 



(20) 



where low resolution bin /3 incorporates the high resolution 
bins a 6 ft, and we set 



Tp a At a 



f x f* R (\,z)d\ 



(21) 



Depending on the galaxy, the final solution obtained with 
the sequence detailed in Section 12.3.21 can be described in 
terms of low-resolution age bins. In this case we should in- 
terpret the recovered mass as the total mass formed during 
the period implied by the width of the bin, but we can- 
not make any conclusions as to when in the bin the mass 
was formed. Similarly, the recovered metallicity for the bin 
should be interpreted as a mass- weighted metallicity for the 
total mass formed in the bin. 



For each galaxy we recover N star formation masses, with 
an associated metallicity, where N is the total number of 
bins, and a maximum of two dust parameters. 



2.4 The models 

The backbone to our g rid of models is the BC 03 set of 
synthetic SSP models (|Bruzual fc Chariot! 120031 ). with a 
Chabrier initial mass func tion (IChabrierll2003l) and Padova 
1994 evolutionary tracks ( Alongi et al.l 19931; iBressan et al.l 



ll993l ; lFagotto et al.lll994al ibl; lGirardi et alJll996h ■ Although 
any set of stellar population models can be used, these 
provide a detailed spectral evolution of stellar populations 
over a suitable range of wavelength, ages and metallicities: 
S(X,t,Z). The models have been normalised to one solar 
mass at the age t = 0. 



2.5 Errors 

The quality of our fits and of our solutions is affected by 
the noise in the data, the noise in the models, and the 
parametrization we choose (which does not reflect the 
complete physical scenario within a galaxy). We aim to 
apply VESPA firstly to SDSS galaxies, which typically have 
a SNR « 20 per resolution element of 3 A, which puts us in 
a regime where the main limitations come from the noise in 
the data. 

To estimate how much noise affects our recovered solutions 
we take a rather empirical approach. For each recovered so- 
lution we create n error random noisy realisations and we 
apply VESPA to each of these spectra. We re-bin each re- 
covered solution in the parametrization of the solution we 
want to analyse and estimate the covariance matrices 
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Figure 3. The evolution of the fit, as VESPA searches for a solution. Sequence should be read from top to bottom. Each line shows 
a stage in the sequence: the left panel shows the input star formation history in the dashed line (red on the online version), and the 
recovered mass fractions on the solid line (black on the online version) for a given parametrization ; the middle panel shows the input 
metallicities in the dashed line (red on the online version), and the recovered metallicities on the solid line (black on the online version); 
the right panel shows the absolute value of the perturbed singular values |u K ■ F| (stars and solid line) and the estimated noise level 
(F) /SNR. In this panel we also show the value of K cr it and the number of non-zero elements of cls in each iteration. The chosen solution 
is the fifth from the top, and indicated accordingly. This galaxy consists of 2 burst events of equal star formation rate - a very young 
and an old burst. It was modelled with a resolution of 3A and a signal-to-noise ratio per pixel of 50. We see the recovery is good but not 
perfect - there is a 1 per cent leakage from the older population - but better than the following solutions, where this bin is split. See text 
in Section l2.3.2l for more details. 



(22) 2.6 Timings 



C{x) a g = ({x a - X a )(xp - Xp)) 

C(Z) afj = ((Z a -Z a )(Z -Z )). (23) 

All the plots in Sections [3] and [4] show error bars derived 
from Cl/a, although it is worth keeping in mind that these 
are typically highly correlated. 



A basic run of VESPA (which consists of roughly 5 runs 
down the sequence detailed in Section 12.3.21 one for each 
value of dust extinction) takes about 5 seconds. If accurate 
error estimations are needed per galaxy, this will add an- 
other one or two minutes to the timing, depending on how 
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accurately one would like to estimate the covariance ma- 
trices, and depending on the number of data points. With 
rierror = 10, a typical SDSS galaxy takes around one minute 
to analyse. 



3 TESTS ON SIMULATED DATA 

We tested VESPA on a variety of synthetic spectra, in 
order to understand its capabilities and limitations. In 
particular, we tried to understand the effect of three factors 
in the quality of our solutions: the input star formation 
history, the noise in the data, and the wavelength coverage 
of the spectrum. We have also looked at the effects of dust 
extinction. Throughout we have modelled our galaxies in a 
resolution of 3A. 

Even though we are aware that showing individual examples 
of VESPA's results from synthetic spectra can be extraor- 
dinarily unrepresentative, we feel obliged to show a few for 
illustration purposes. We will show a typical result for most 
of the cases we present, but we also define some measure- 
ments of success, so that the overall performance of VESPA 
can be tracked as we vary any factors. We define 



and 



Gz = 



(24) 



(25) 



where x 1 ^ and Z a are the total mass and correspondent 
metallicity in bin a (re-binned to match the solution's 
parametrization if necessary), and cu a is the flux contribu- 
tion of population of age t a . G x and Gz are a flux- weighted 
average of the total absolute fractional errors in the solu- 
tion, and give an indication of how well VESPA recovers 
the most significant parameters. A perfect solution gives 
G x = Gz = 0. It is also worth noting that this statistic 
does not take into account the associated error with each 
recovered parameter - deviations from the true solution are 
usually expected given the estimated covariance matrices. 
We will also show how these factors affect the recovered to- 
tal mass for a galaxy. In all cases we have re-normalised the 
total masses such that total input mass for each galaxy is 1. 

3.1 Star formation histories 



We present here some results for synthetic spectra with 
two different star formation histories. All of the spectra 
in this section were synthesised with a SNR per pixel of 
50, and we initially fit the very wide wavelength range 
A e [1000, 9500]A. 

We choose two very difference cases: firstly a star formation 
history of dual bursts, with a large random variety of 
burst age separations and metallicities (where we set the 
star formation rate to be 10 solar masses per Gyr in all 
bursts). Secondly, we chose a SFH with an exponentially 
decaying star formation rate: SFR oc exp(7t a ), where t a is 



0.3 - 



0.' - 



5 10 15 20 

Number of recovered non-zero parameters 

Figure 8. The recovered number of non-zero parameters in 50 
galaxies with an exponentially decaying star formation history, 
using different wavelength ranges: A e [1000, 9500] A(solid line) 
and A G [3200, 9500] A(dashed line). Please note that these corre- 
spond to the total number of non-zero components in the solution 
vector c K and not to the number of recovered stellar populations. 



the age of the bin in lookback time in Gyr. Here we show 
results for 7 = 0.3 Gyr -1 . Rather than being physically 
motivated, our choice of 7 reflects a SFH which is not 
too steep as to essentially mimic a single old burst, but 
which is also not completely dominated by recent star 
formation. In all cases the metallicity in each bin is ran- 
domly set. Figure |3] shows a typical example from each type. 

Figure [5] shows the distribution of G x , Gz and of the recov- 
ered total masses for a sample of 50 galaxies. We see differ- 
ences between the two cases. Firstly, in dual bursts galaxies, 
we seem to do better in recovering data from significant in- 
dividual bins, but worse in overall mass. This reflects the 
fact that G x is dominated by the fractional errors in the 
most significant bins, but the total mass can be affected by 
small flux contributions in old bins which can have large 
masses. On the other hand, with an exponentially decaying 
star formation rate, we do worse overall (although this is 
mainly a reflection that more bins have significant contribu- 
tions to the flux) but we recover the total mass of the galaxy 
exceptionally well. 



3.2 Wavelength range 

Wavelength range is an important factor in this sort of 
analysis, as different parts of the spectrum will help to break 
different degeneracies. Since we are primarily interested in 
SDSS galaxies, we have studied how well VESPA does in 
the more realistic wavelength range of A € [3200, 9500] A. 

Figure [6] shows the results for the same galaxies shown in 
Figure [4] obtained with the new wavelength range. In these 
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Figure 4. Two examples of VESPA's analysis on synthetic galaxies. The top panels show the original spectrum in the dark line (black 
in the online version) and fitted spectrum in the lighter line (red in the online version ). The middle panels show the input (dashed, red) 
and the recovered (solid, black) star formation rates and the bottom panel shows the input (dashed, red) and recovered (solid, black) 
metallicities per bin. Note that even though many of the recovered metallicities are wrong, these tend to correspond to bins with very 
little star formation, and are therefore virtually unconstrained. 




Figure 5. The distribution of G x , Gz and total mass recovered for 50 galaxies with a SNR per pixel of 50. Solid lines correspond to 
dual burst and dashed lines to exponentially decaying ones. See text in Section [3T] for details. 



particular cases, we notice a more pronounced difference in 
the dual bursts galaxy, but looking at a more substantial 
sample of galaxies shows that this is not generally the case. 
Figure [7] shows G x , Gz and total mass recovered for 50 
exponentially decaying star formation history galaxies, with 
a signal to noise ratio of 50 and the two different wavelength 
ranges. We do not see a largely significant change in both 
cases, and we observe a less significant difference in the 
dual bursts galaxies (not plotted here). 

We find it instructive to keep track of how many parame- 
ters we recover in total, as we change any factors. Figure 
[8] shows an histogram of the total number of non-zero 



parameters we recovered from our sample galaxies with 
exponentially-decaying star formation histories and both 
wavelength ranges. Note that these are the components of 
the solution vector c K which are non-zero - they do not 
represent a number of recovered stellar populations. In this 
case there is a clear decrease in the number of recovered 
parameters, suggesting a wider wavelength range is a useful 
way to increase resolution in parameter space. 
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Figure 6. Same galaxies as in Figure [4] but results obtained by using a smaller wavelength range. The goodness-of-fit in data space is 
still excellent, but it becomes more difficult to break certain degeneracies. 
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Figure 7. The distribution of G x , Gz and total mass recovered for 50 galaxies with a SNR per pixel of 50 and two different wavelength 
coverage. Solid line corresponds to A G [1000, 9500] A and dashed line to A £ [3200, 9500] A. 



3.3 Noise 

It is of interest to vary the signal-to-noise ratio in the 
synthetic spectra. We have repeated the studies detailed 
in the two previous sections with varying values of noise, 
and we investigate how this affects both the quality of the 
solutions and their resolution in parameter space. 

Figure [10] shows how the recovered number of parameters 
changed by increasing the noise in the galaxies with an 
exponentially decaying star formation rate and wide wave- 
length range. In this case the increase in the noise leads to a 
significant reduction of the number of parameters recovered 
for each galaxy. This behaviour is equally clear for different 
star formation histories and different wavelength coverage, 
and is directly caused by the stopping criterion defined in 



Section 13.31 

The quality of the solutions is also affected by this increase 
in noise, as can be seen in figure|9] where we have plotted G x , 
Gz and the total recovered mass for two different values of 
SNR. The quality of the solutions decreases with the higher 
noise levels, as is to be expected. However, a more interest- 
ing question to ask is whether this decrease in the quality 
of the solutions would indeed be more pronounced without 
the SVD stopping criterion. Figure [TT] shows a comparison 
between G x obtained as we have described and obtained 
without any stopping mechanism (so letting our search go 
to the highest possible resolution and taking the final so- 
lution) for 50 galaxies with an exponentially decaying star 
formation history and a signal-to-noise ratio of 20. The re- 
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Figure 9. The distribution of G x , Gz and total mass recovered for 50 galaxies with an exponential decaying star formation history and 
different signal-to-noise ratios. Solid lines correspond to SNR = 50 and dashed lines to SNR =20. See text in Section 13.31 for details. 
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Figure 10. The recovered number of non-zero parameters as we 
change the noise in the data from 50 (solid line) to 20 (dashed 
line), in a sample of galaxies with an exponentially decaying star 
formation rate. Please note that these correspond to the total 
number of non-zero components in the solution vector c K and 
not to the number of recovered stellar populations. 



Figure 11. Testing the SVD stopping criterion. Plots show good- 
ness of fit G x for the solution of 50 galaxies obtained with and 
without the SVD stopping criterion. We see that by recovering 
only as many parameters as the data warrants gives improved 
parameter estimation in almost all cases, and a striking improve- 
ment in many. 



suits show clearly that there is a significant advantage in 
using the SVD stopping criterion. Naturally, the goodness 
of fit in data space is consistently better as we increase the 
number of parameters but this improvement is illusory - the 
parameter recovery is worse. This is exactly the expected be- 
haviour - we choose to sacrifice resolution in parameter space 
in favour of a more robust solution - even though naively one 
could think a lower x 2 solution would indicate a better solu- 
tion. The significance of this improvement changes with the 
amount of noise and wavelength range of the data (and to 
a lesser extent with type of star formation history) but we 
observed an improvement in all cases we have studied. 



As expected, further decreasing the signal-to-noise ratio 
leads to a further degradation of the recovered solutions. 
This is accompanied by a suitable increase in the error bars 
and correlation matrices, but in cases of a SNR~ 10 and less 
it becomes very difficult to recover any meaningful informa- 
tion from individual spectra. 

3.4 Dust 

In this section we use simulated galaxies to study the effect 
of dust in our solutions. As explained in Section 12.1.11 
due to the non-linear nature of the problem, we cannot 
include dust as one of the free parameters analysed by 
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Figure 12. Testing the recovery of T-L for 50 galaxies with 
a exponentially-decaying star formation history (triangles) and 
50 galaxies formed with a random combination of dual bursts 
(stars). The input values are randomly chosen and continuously 
distributed between and 2. The recovered values are chosen from 
a tabulated grid between and 4. 



when using the more sophisticated model, we noted that the 
mean error on r v SM on a subsample of dual-burst galaxies 
(synthesised as explained in section I3TT1 but chosen to have 
young star formation) was reduced from 35 to 28 per cent. 
This simple test also revealed that we are less likely to 
underestimate the mass of young populations by allowing 
an extra dust component, but that we are also introducing 
an extra degeneracy, especially so in the case of faint young 
populations. However, we feel that the two-parameter dust 
model brings more advantages than disadvantages, with 
the caveat being that dusty young populations can be 
poorly constrained. In any case, we note that each galaxy 
is always analysed with a one-parameter model before 
being potentially analysed with a two-parameter model, 
and both solutions are kept and always available for analysis. 

Finally, our test also partly justifies our choice to first 
run a single dust component model and only apply a two- 
component model if we detect stars in the first two bins - 
we find that although a one-component model might un- 
derestimate the amount of young stars, it does not fail 
to detect them. We repeated a similar test in real data, 
by analysing the same sample with and one- and a two- 
parameter dust model. We found similar results, with a one- 
parameter model failing to yield star formation in young bins 
only around 1 per-cent of the time (compared to the two- 
parameter model), and only in cases where the contribution 
of the light from the young populations was very small (of 
the order of 1 to 2 per cent). 



SVD. Instead, we fit for a maximum of two dust parameters 
using a brute-force approach which aims to minimise % 2 m 
data-space by trying out a series of values for t v sm and Ty c . 

For each galaxy we assign random values of t v sm £ [0, 2] 
and Ty C 6 [1,2] and we are interested in how well we 
recover these parameters and any possible degeneracies. 

Figure [12] shows the input and recovered values for t v sm for 
galaxies with a signal-to-noise ratio of 50, and which were 
analysed using the wavelength range A G [3200, 9500]A. We 
show results for two different cases of star formation his- 
tory: 50 galaxies with an exponentially-decaying SFR and 
50 galaxies formed by dual-bursts. We observe a good recov- 
ery of t v sm in both cases, especially at low optical depths. 

However, we mostly observe a poor recovery of Ty C , 
especially at high optical depths. This is unsurprisingly 
flagging up a certain level of degeneracy between mass and 
degree of extinction, which gets worse as the optical depth 
increases. Essentially, it becomes difficult to distinguish 
between a highly obscured massive population and a less 
massive population surrounded by less dust. It is worth 
keeping in mind that young populations are affected by 
both dust components simultaneously, and generally, even 
though the recovery of the second dust parameter may 
not be accurate, it allows for a better estimation of the 
dominant dust component. 

This can be tested by simulating galaxies on a two- 
component dust model and by analysing them using both a 
single component model, and a two-component model. E.g., 



4 RESULTS 

In this section we present some results obtained by applying 
VESPA to galaxies in the SDSS. Our aim is to analyse these 
galaxies, and to produce and publish a catalogue of robust 
star formation histories, from which a wealth of information 
can then be derived. We leave this for another publication, 
but we present here results from a sub-sample of galaxies, 
which we used to test VESPA in a variety of ways. 

4.1 Handling SDSS data 

Prior to any analysis, we processed the SDSS spectroscopic 
data, so as to accomplish the desired spectral resolution 
and mask out any wanted signal. 

The SDSS data-files supply a mask vector, which flags 
any potential problems with the measured signal on a 
pixel-by-pixel basis. We use this mask to remove any 
unwanted regions and emission lines. In practical terms, we 
ignore any pixel for which the provided mask value is not 
zero. 

The BC03 synthetic models produce outputs at a resolution 
of 3A, which we convolve with a Gaussian velocity disper- 
sion curve with a stellar velocity ay = 170kms _1 , this being 
a typical value for SDSS galaxies. We take the models' 
tabulated wavelength values as a fixed grid and re-bin the 
SDSS data into this grid, using an inverse- variance weighted 
average. We compute the new error vector accordingly. Note 
that the number of SDSS data points averaged into any new 
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bin is not constant, and that the re-binning process is done 
after we have masked out any unwanted pixels. Additionally 
to the lines yielded by the mask vector, we mask out the 
following emission line regions in every spectrum's rest- 
frame wavelength range: [5885-5900, 6702-6732, 6716-6746, 
6548-6578, 6535-6565, 6569-6599, 4944-4974, 4992-5022, 
4846-4876, 4325-4355, 4087-4117, 3711-3741, 7800-11000] A. 

These re-binned data- and noise-vectors are essentially the 
ones we use in our analysis. However, since the linear algebra 
assumes white-noise, we pre- whiten the data and construct 
a new flux vector Fj — Fj/oj, which has unit variance, 
a'j — l,Vj, and a new model matrix = Aijjoj. 



4.2 Duplicate galaxies 

There are a number of galaxies in the SDSS database 
which have been observed more than once, for a variety 
of reasons. This provides an opportunity to check how 
variations in observation-dependent corrections affect the 
results obtained by VESPA. 

We have used a subset of the sample of duplicate objects in 
iBrinchmann et all ([2004 13 to create two sets of oservations 
for 2000 galaxies, which we named list A and list B. 
We are interested in seeing how the errors we estimate 
for our results compare to errors introduced by intrinsic 
variations caused by changing the observation conditions 
(such as quality of the spectra, placement of the fibre, sky 
subtraction or spectrophometric calibrations). 

Figure [13] shows the average star formation fraction as a 
function of lookback time for both sets of observation. The 
error bars showed are errors on the mean. We see no signs 
of being dominated by systematics when estimating the 
star formation fraction of a sample of galaxies. 



Figure [14] shows the total stellar mass obtained for a set of 
500 galaxies in both observations (details of how we estimate 
the total stellar mass of a galaxy are included in section |4~41 . 
The error bars are obtained directly from the estimated 
covariance matrix C(x) (equation I22[) . Even though most of 
the galaxy duplicates produce mass estimates in agreement 
with each other given the error estimates, a minority 
does not. Upon inspection, these galaxies show significant 
differences in their continuum, but after further investiga- 
tion it remains unclear what motivates such a difference. 
The simplest explanation is that the spectrophotometric 
calibration differs significantly between both observations, 
and that might have been the reason the plate or object 
was re-observed. Whatever the reason however, the clear 
conclusion is that stellar mass estimates are highly sensitive 
to changes in the spectrum continuum, and the errors we es- 
timate from the covariance matrix alone might be too small. 



We did not find any signs of a systematic bias in any of the 
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Figure 13. Average star formation fraction as a function of look- 
back time for the 2000 galaxies in list A (solid line) and list B 
(dashed line). The error bars shown are the errors bars on the 
mean for each age bin. We show only the errors from list A to 
avoid cluttering - the errors from list B are of similar amplitude. 
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Figure 14. Total stellar mass recovered for two sets of observa- 
tions of 500 galaxies in the main galaxy sample. The error bars 
are calculated from C(x). 



analysis we carried out. 



4.3 Real fits 

In this section we discuss the quality of the fits to SDSS 
galaxies obtained with VESPA. 



1 Available at http://www.mpa-garching.mpg.de/SDSS/ 



As explained in Section 2, VESPA finds the best fit solution 
in a x 2 sense for a given parametrisation, which is self- 
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Figure 15. The distribution of reduced values of \ 2 f° r a sample 
of 360 galaxies analysed by VESPA. 

regulated in order to not allow an excessive number of fitting 
parameters. We have shown that this self-regularization 
gives a better solution in parameter space (Figure Hip , 
despite often not allowing the parametrization which would 
yield the best fit in data space (Figure [3}. However, our 
aim is still to find a solution which gives a good fit to the 
real spectrum. Figure \T5\ shows the 1-point distribution of 
reduced values of reduced \ 2 f° r 1 plate of galaxies. This 
distribution peaks at around Xr educed = 1-3, and figure [TrJl 
shows a fit to one of the galaxies with a typical value of 
goodness of fit. 

It is worth noting that the majority of the fits which are 
most pleasing to the eye, correspond to the ones with a 
high signal to noise ratio and high value of reduced \ 2 . One 
would expect the best fits to come from the galaxies with 
the best signal. However, we believe the fact that they do 
not is not a limitation of the method, but a limitation of 
the modelling. There are a number of reasons why VESPA 
would be unable to produce very good fits to the SDSS 
data. One is the adoption of a single velocity dispersion 
(170 kms -1 ) which could easily be improved upon at the 
expense of CPU time. However, the dominant reason is 
likely to be lack of accuracy in stellar and dust modelling 
- whereas BC03 models can and do reproduce a lot of the 
observed features, it is also well known that this sucess 
is limited as there are certain spectral features not yet 
accurately modelled, or even modelled at all. There are 
similar deficiencies in dust models and dust extinction 
curves. The effect of the choice of modelling should not be 
overl ooked, and we refer t he reader to a discussion in Section 
4.5 of lPanter et"afl ((2006)), where these issues are discussed. 



4.4 VESPA and MOPED 

In this Section we take the opportunity to compare the 
results from VESPA and MOPED, obtained from the same 



Figure 17. The recovered average star formation history for the 
821 galaxies as recovered by VESPA (solid line) and MOPED 
(dashed line). Both were initially normalised such that the sum 
over all bins is 1, and the MOPED line was then adjusted by 
11/16 to account for the different number of bins used in each 
method, to facilitate direct comparison. 



sample of galaxies. The VESPA solutions used here are 
obtained with a one-parameter dust model, to allow a more 
fair comparison between the two methods. Both methods 
make similar assumptio ns regarding stellar models, but 
MOPED uses an LMC l|Gordon et alj |2003( ) dust extinc- 
tion curve, and single screen modelling for all optical depths. 

Our sample consist o f two plates from the SDSS DR3 
|Abazaiian et alj|2005h (plates 0288 and 0444), from which 
we analyse a total of 821 galaxies. We are mainly inter- 
ested in comparing the results in a global sense. MOPED in 
its standard configuration attempts to recover 23 parame- 
ters (11 star formation fractions, 11 metallicities and 1 dust 
parameter), so we might expect considerable degeneracies. 
Indeed, in the past the authors of MOPED have cautioned 
against using it to interpret individual galaxy spectra too 
precisely. We have observed degeneracies between adjacent 
bins in MOPED, but on the other hand a typical MOPED 
solution has many star formation fractions which are es- 
sentially zero, so the number of significant contributions is 
always much less than 23. 

Figure [17] shows the recovered average star formation 
history for the 821 galaxies using both methods. In the 
case of VESPA, solutions parametrized by low-resolution 
bins had to be re-parametrized in high-resolution bins, so 
that a common grid across all galaxies could be used. This 
was done using the weights given by (|2ip . The lines show a 
remarkably good agreement between the two methods. 

Having recovered a star formation history for each galaxy, 
one can then estimate the stellar mass of a galaxy. We calcu- 
lated this quantity for all galaxies using the solutions from 
both methods, and with similar assumptions regarding cos- 
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Figure 16. Typical fit to a galaxy from the SDSS. The dark line is the real data (arbitrary normalisation), and the lighter line (red on 
the online version) is VESPA's fit to the data. 



mological parameters and fibre-size corrections. Explicitly, 
we have done the following: 

(i) We converted from flux to lumi nosity assuming the set 
of cosmological parameters given bv lSpergel et alj ((2003). 

(ii) We recovered the initial mass in each age bin using 
each method. 

(iii) We calculated the remaining present-day mass for 
each population after recycling processes. This information 
is supplied by the synthetic stellar models, as a function of 
age and metallicity. 

(iv) We summed this across all bins to calculate the total 
stellar present-day mass in the fibre aperture, M. 

(v) We corrected for the aperture size by scaling up the 
mass to Mstellar using the petrosian and fibre magnitudes 
in the z-band, M p (z) and Mf(z), with: M ste uar = M x 

IQ0 A[M p (z)-M f (z)} 

Figure [18] shows the recovered galaxy masses as recovered 



from MOPED and from VESPA. We see considerable agree- 
ment between VESPA and MOPED. Over 75 per cent of 
galaxies have 0.5 < Mvespa/ Mmoped < 1-5. There is a 
tail of around 10 per cent of galaxies where VESPA recov- 
ers 2 to 4 times the mass recovered by MOPED. The main 
reason for this difference is in the dust model used - we find 
a correlation between dust extinction and the ratio of the 
two mass estimates. This again reflects the fact that total 
stellar mass estimates are highly sensitive to changes in the 
spectrum continuum (see also section l4~2jl . 

Our sub-sample of 821 includes galaxies with a wide range 
of signal-to-noise ratios, star formation histories and even 
wavelength range (mainly due to each galaxy having dif- 
ferent masks applied to it, according to the quality of the 
spectroscopic data). Figure [T^l shows the number of recov- 
ered non-zero parameters in the sample, using VESPA. As 
an average, it falls below the synthetic examples studied in 
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Figure 18. Galaxy stellar mass (in units of solar masses) as re- 
covered by VESPA and MOPED for a sub-sample of 821 SDSS 
galaxies. The small percentage of galaxies with significantly larger 
VESPA masses have large extinction. The difference is accounted 
for by the fact that MOPED and VESPA use different dust mod- 
els. 



Figure 19. Number of non-zero parameters in solutions recov- 
ered from 821 SDSS galaxies with VESPA. Please note that these 
correspond to the total number of non-zero components in the 
solution vector c K and not to the number of recovered stellar 
populations. For a information about the number of recovered 
populations see Figure 14. 



Section 3. This is not surprising, though, as each galaxy will 
have an unique and somewhat random combination of char- 
acteristics which will lead to a different number of param- 
eters being recovered. The total combination of these sets 
of characteristics would be impossible to investigate using 
the empirical method described in Section 3, and here lies 
the advantage of VESPA of dynamically adapting to each 
individual case. Also important to note is the fact that the 
wavelength coverage is normally not continuous in an SDSS 
galaxy, due to masked regions. This was not modelled in 
Section 3, and is likely to further reduce the number of re- 
covered parameters in any given case. 

Perhaps more useful is to translate this number into a 
number of recovered significant stellar populations for each 
galaxy. We define a significant component as a stellar pop- 
ulation which contributes 5 per cent or more to the total 
flux. Figure [201 shows the distribution of the number of sig- 
nificant components for our sub-sample of galaxies, as re- 
covered by MOPED and VESPA. ft is interesting to note 
that both methods recover on average a similar amount 
of components, even though MOPED has no explicit self- 
regularization mechanism, as VESPA clearly does. 



5 CONCLUSIONS 

We have developed a new method to recover star formation 
and metallicity histories from integrated galactic spectra 
- VESPA. Motivated by the current limitations of other 
methods which aim to do the same, our goal was to develop 
an algorithm which is robust on a galaxy-by-galaxy basis. 
VESPA works with a dynamic parametrization of the 
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Figure 20. The distribution of the total number of recovered 
stellar populations which contribute 5 per-cent or more to the 
total flux of the galaxy, as recovered from MOPED (dashed line) 
and VESPA (solid line). 
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star formation history, and is able to adapt the number 
of parameters it attempts to recover from a given galaxy 
according to its spectrum. In this paper we tested VESPA 
against a series of idealised synthetic situations, and against 
SDSS data by comparing our results with those obtained 
with the well-established code, MOPED. 

Using synthetic data we found the quality and resolution 
of the recovered solutions varied with factors such as type 
of star formation history, noise in the data and wavelength 
coverage. In the vast majority of cases, and within the 
estimated errors and bin-correlations, we observed a reliable 
reproduction of the input parameters. As the signal-to-noise 
decreases, it becomes increasingly difficult to recover robust 
solutions. Whereas our method cannot guarantee a perfect 
solution, we have shown that the self-regularization we im- 
posed helped obtain a cleaner solution in an overwhelming 
majority of the cases studied. 

On the real data analysis, we have studied possible effects 
from systematics using duplicate observations of the same 
set of galaxies, and have also compared VESPA's to 
MOPED 's results obtained using the same data sample. We 
found that in the majority of cases our results are robust 
to possible systematics effects, but that in certain cases 
and particularly when calculating stellar masses, VESPA 
might underestimate the mass errors. However, we found 
no systematic bias in any of our tests. We have also shown 
that VESPA's results are in good agreement with those 
of MOPED for the same sample of galaxies. VESPA and 
MOPED are two fundamentally different approaches to 
the same problem, and we found good agreement both in 
a global sense by looking at the average star formation 
history of the sample, and in an individual basis by looking 
at the recovered stellar masses of each galaxy. VESPA 
typically recovered between 2 to 5 stellar populations from 
the SDSS sample. 

VESPA's ability to adapt dynamically to each galaxy and 
to extract only as much information as the data warrant is 
a completely new way to tackle the problem of extracting 
information from galactic spectra. Our claim is that, for 
the most part, VESPA's results are robust for any given 
galaxy, but our claim comes with two words of caution. The 
first one concerns very noisy galaxies - in extreme cases 
(SNR«10 or less, at a resolution of 3A), it becomes very 
difficult to extract any meaningful information from the 
data. This uncertainty is evident in the large error bars 
and bin-correlations, and the solutions can be essentially 
unconstrained even at low-resolutions. We are therefore 
limited when it comes to analysing individual high-noise 
galaxies, which is the case of many SDSS objects. Our 
second word of caution concerns the stellar models used to 
analyse real galaxies - any method can only do as well as the 
models it bases itself upon. We are limited in our knowledge 
and ability to reproduce realistic synthetic models of stellar 
populations, and this is inevitably reflected in the solutions 
we obtain by using them. On the plus side, VESPA works 
with any set of synthetic models and can take advantage of 
improved versions as they are developed. 

VESPA is fast enough to use on large spectroscopic sam- 



ples (a typical SDSS galaxy takes 1 minute on an average 
workstation), and we are in the process of analysing SDSS's 
Data Release 5 (DR5), which consists of roughly half a mil- 
lion galaxies. Our first aim is to publish and exploit a cata- 
logue of robust star formation histories, which we hope will 
be a valuable resource to help constrain models of galaxy 
formation and evolution. 
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